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Abstract 

Mapping expression Quantitative Trait Loci (eQTLs) represents a powerful and widely-adopted approach 

to identifying putative regulatory variants and linking them to specific genes. Up to now eQTL studies 
have been conducted in a relatively narrow range of tissues or cell types. However, understanding the 
biology of organismal phenotypes will involve understanding regulation in multiple tissues, and ongoing 
studies are collecting eQTL data in dozens of cell types. Here we present a statistical framework for 
powerfully detecting eQTLs in multiple tissues or cell types (or, more generally, multiple subgroups). 
The framework explicitly models the potential for each eQTL to be active in some tissues and inactive in 
others. By modeling the sharing of active eQTLs among tissues this framework increases power to detect 
eQTLs that are present in more than one tissue compared with "tissue-by-tissue" analyses that examine 
each tissue separately. Conversely, by modeling the inactivity of eQTLs in some tissues, the framework 
allows the proportion of cQTLs shared across different tissues to be formally estimated as parameters of a 
model, addressing the difficulties of accounting for incomplete power when comparing overlaps of cQTLs 
identified by tissue-by-tissuc analyses. Applying our framework to re-analyze data from transformed B 
cells, T cells and fibroblasts we find that it substantially increases power compared with tissuc-by-tissue 
analysis, identifying 63% more genes with eQTLs (at FDR=0.05). Further the results suggest that, in 
contrast to previous analyses of the same data, the majority of eQTLs detectable in these data are shared 
among all three tissues. 

Author Summary 

Genetic variants that are associated with gene expression arc known as expression Quantitative Trait Loci, 
or eQTLs. Many studies have been conducted to identify eQTLs, and they have proven an effective tool 
for identifying putative regulatory variants and linking them to specific genes. Up to now most studies 
have been conducted in a single tissue or cell-type, but moving forward this is changing, and ongoing 
studies are collecting data aimed at mapping eQTLs in dozens of tissues. Current statistical methods 
are not able to fully exploit the richness of these kinds of data, taking account of both the sharing and 
differences in eQTLs among tissues. In this paper we develop a statistical framework to address this 
problem, to improve power to detect eQTLs when they are shared among multiple tissues, and to allow 
for differences among tissues to be estimated. Applying these methods to data from three tissues suggests 
that sharing of cQTLs among tissues may be substantially more common than it appeared in previous 
analyses of the same data. 
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Introduction 

Regulatory variation plays an essential role in the genetics of disease and other phenotypes as well as 
in evolutionary change [l}|3]. However, in sharp contrast to nonsynonymous variants in the human 
genome, which can now be identified with great accuracy, it remains extremely difficult to know which 
variants in the genome may impact gene regulation in any given tissue or cell type. [Henceforth we use 
"tissue" for brevity, but everything applies equally to cell types.] Expression QTL mapping (e.g. |4||6] 
represents a powerful approach for bridging this gap, by allowing regulatory variants to be identified, 
and linked to specific genes. Indeed, numerous studies (e.g., [7j|8]) have shown highly significant overlaps 
between eQTLs and SNPs associated with organismal-level phenotypes in genome- wide association studies 
(GWAS), suggesting that a large fraction of GWAS associations may be due to variants that affect gene 
expression. 

Ultimately, understanding the biology of organismal phenotypes, such as diseases, is likely to require 
understanding regulatory variation in many different tissues ( [9{ |10| ). For example, if regulatory variants 
differ across tissues, then, in understanding GWAS hits, and using them to understand the biology of 
disease, we would like to know which variants are affecting which tissues. At a more fundamental level, 
identifying differential genetic regulation in different tissues could yield insights into the basic biological 
processes that influence tissue differentiation. To date, eQTL studies have been performed in a relatively 
narrow range of tissue types. However, this is changing quickly: for example, the NIH "Genotype- 
Tissue Expression" (GTEx) project aims to collect expression and genotype data in 30 tissues across 
900 individuals. Motivated by this, here we describe and illustrate a statistical framework for mapping 
eQTLs in expression data on multiple tissues. 

While statistical methods for identifying eQTLs in a single tissue or cell type are now relatively mature 
(e.g. [IT]) current analytic tools are limited in their ability to fully exploit the richness of data across 
multiple tissues. In particular, available methods fall short in their ability to jointly analyze data on all 
tissues to maximize power, while simultaneously allowing for differences among eQTLs present in each 



tissue. Indeed relatively few papers have considered the problem. The simplest approach (e.g. 12 13]) 
is to analyze data on each tissue separately ("tissue-by-tissue" analysis), and then to examine overlap of 
results among tissues. However, this fails to leverage commonalities among tissues to improve power to 
detect shared eQTLs. Furthermore, although examining overlap of eQTLs among tissues may appear a 
natural approach to examining heterogeneity, in practice interpretation of results is complicated by the 
difficulty of accounting for incomplete power. Both |14 and |13 provide approaches to address this, but 
only for pairwise comparisons of tissues. 

Compared with tissue-by-tissue analysis, joint analysis of multiple tissues has the potential to increase 



power to identify eQTLs that have similar effects across tissues. Both 15 and 16 conduct such joint 
analyses - the first using ANOVA, and the second using a weighted Z-score meta-analysis - and 16 
confirm that their joint analysis has greater power than tissue-by-tissue analysis. The ANOVA and 
Z-score methods each have different advantages. The ANOVA framework has the advantage that, by 
including interaction terms, it can be used to investigate heterogeneity in eQTL effects among tissues. 
Gerrits et al. ( [l5]) use this to identify eQTLs that show significant heterogeneity, and then classify these 
eQTLs, post-hoc, into different types based on estimated effect sizes. The weighted Z-score method has 
the advantage that, unlike ANOVA, it allows for different variances of expression levels in different 
tissues (which are likely to occur in practice). However, it does not so easily allow for investigation 
of heterogeneity; Fu et al. ( [TBj) assess heterogeneity for pairs of tissues by using a resampling-based 
procedure to assess the significance of observed differences in Z scores. 

Here we introduce a statistical framework for the joint analysis of eQTLs among multiple tissue types, 
that combines advantages of some of the methods above, as well as introducing some new ones. In brief, 
our framework integrates recently-developed GWAS meta-analysis methods that allow for heterogeneity 



of effects among groups 17 -20 , into a hierarchical model (e.g. |21[|22|) that combines information across 



genes to estimate the relative frequency of patterns of eQTL sharing among tissues. Like ANOVA, our 
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approach allows investigation of heterogeneity among several tissues, not just pairs of tissues. However, 
in contrast to ANOVA, our framework allows for different variances in different tissues. Moreover, unlike 
any of the methods described above, our framework explicitly models the fact that some tissues may share 
eQTLs more than others, and estimates these patterns of sharing from the data. Our methods also allow 
for intra-individual correlations when samples are obtained from a common set of individuals. While we 
focus here on comparing and combining information across different tissue types, our framework could 
be applied equally to comparing and combining across other units, e.g. different experimental platforms, 
multiple datasets on the same tissue types, or data on individuals from different populations. 

The remainder of the paper is as follows. After providing a brief overview of our framework, we use 
simulations to illustrate its power compared to other methods, and then apply it to map eQTLs, and 
assess heterogeneity among tissues, using data from Fibroblasts, LCLs and T-cells ( [12|). Consistent 
with results from 16 , we show that our joint analysis framework provides a large gain in power compared 
with a tissue-by-tissue analysis. Furthermore, compared with previous analyses of these data, we find a 
much higher rate of tissue-consistent eQTLs. 



Results 

Methods Overview 

Consider mapping eQTLs in S tissues. In our applications here the expression data are from micro-arrays, 
and so we assume a normal model for the expression levels, suitably-transformed. (These methods can 
also be applied to RNA-seq data after suitable transformation; see Discussion). That is, in each tissue, 
s = 1, . . . , 5, we model the potential association between a candidate SNP and a target gene by a linear 
regression: 

Vsi = fJ-s+ Ps9i + f-si with esi ^ 7V(0, cr^), (1) 

where ysi denotes the observed expression level of the target gene in tissue s for the i*'* individual, /i^ the 
mean expression level of this gene in tissue s, j3s the effect of a candidate SNP on this gene expression in 
tissue s, Qi the genotype of the i*'* individual at the SNP (coded as 0,1 or 2 copies of a reference allele) 
and esi the residual error for tissue s and individual i. Note that the subscript s on residual variance 
indicates that we allow the residual variance to be different in each tissue. In addition, when tissues are 
sampled from the same set of individuals, we allow that the residual errors eii, . . . , tsi may be correlated 
(with the correlation matrix to be estimated from the data). 

The primary questions of interest are whether the SNP is an eQTL in any tissue, and, if so, in which 
tissues. To address these questions we use the idea of a "configuration" from |19j|20|. A configuration 
is a binary vector 7 = (71, . . . ,75) where 7s G {0, 1} indicates whether the SNP is an eQTL in tissue 
s. If 7s = 1 then we say the eQTL is "active" in tissue s. The "global null hypothesis", Hq^ that the 
SNP is not an eQTL in any tissue, is therefore 7 = (0, ... ,0). Every other possible value of 7 can be 
thought of as representing a particular alternative hypothesis. For example, 7 = (1,...,1) represents 
the alternative hypothesis that the SNP is an eQTL in all S tissues, and 7 = (1, 0, . . . , 0) represents the 
alternative hypothesis that the SNP is an eQTL in just the first tissue. 

Our aim is to perform inference for 7. A natural approach is to specify a probability model, P(data | 7), 
being the probability of obtaining the observed data if the true configuration were 7, and then perform 
likelihood-based inference for 7. The support in the data for each possible value of 7, relative to the null 
Hq, is quantified by the likelihood ratio, or Bayes Factor (BF, [23]): 



P(data I true configuration is 7) 
P(data |i?o) 



^^7 = n/j^^^ I Tj \ ■ (2) 



Specifying these likelihoods requires assumptions about P(/3|7), the distribution of the effect sizes /? for 
each possible configuration 7 (as well as less crucial assumptions about nuisance parameters such as and 
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(Is). Of course, if 7s = then /3s = by definition, but for the tissues where 7s = 1 various assumptions 
are possible - for example, one could assume that the effect /3s is the same in all these tissues, or allow it 
to vary among tissues. Here we use a flexible family of distributions, P(/?|7, 0) (see Methods), where the 
hyper-parameters 6 can be varied to control both the typical effect size, and the heterogeneity of effects 
across tissues (see below). 

The value of BF^ measures the support in the data for one specific alternative configuration 7, com- 
pared against the null hypothesis Hq. To account for the fact that there are many possible alternatives, 
the overall strength of evidence against at the candidate SNP is obtained by "Bayesian Model Av- 
eraging" (BMA), which involves averaging BF^ over the possible alternative configurations 7, weighting 
each by its prior probability, 1]^: 

_ P(data I Hp false) _ ^ ^.^ 
^^^^^ - P(data I Ho true) " ^ ^^''^ 

Further, under an assumption of at most one eQTL per gene, the overall evidence against Hq for the 
entire gene (i.e. that the gene contains no eQTL in any tissue) is given by averaging BFbma across all 



candidate SNPs 24 . In either case, at either the SNP or gene level, large values of BFbma constitute 
strong evidence against Hq. BFbma has a direct Bayesian interpretation as the strength of the evidence 
against Hq, but here we also use it as a frequentist test statistic ( [24][25j), assessing significance by 
permutation or simulation. The latter has the advantage that p-values and q-values obtained in this way 
are "valid" even if not all the prior assumptions are exactly correct. 

Note that BFbma depends on the choice of {9,7]), and the power of BFbma ^ ^ test statistic is 
expected to depend on how well this choice of these hyper-parameters captures the range of alternative 
scenarios present in the data. Here we make use of three different choices: 

A "data-driven" choice, where the hyper-parameters are estimated from the data using a hierarchical 



model (HM, 26 ) that combines information across all genes. We use BFbma to denote this choice 



• A "default" choice, which chooses ry to cover a wide range of different possible alternative configu- 
rations, and 6 is set to allow modest heterogeneity. We use BFbma to denote this choice. 

• A "lite" choice, which puts weight only on the most extreme configurations (where the eQTL 
is active in only one tissue, or in all tissues), but compensates by setting 9 to allow for more 
heterogeneity. We use BFBMAiitc to denote this choice. 

Each of these choices has something to recommend it. The first, being data driven, is the most attractive 
in principle, but also the most complex to implement. The default choice is simpler to implement, and is 
included partly to demonstrate that one does not have to get the hyper-parameter values exactly "right" 
for BFbma to be a powerful test statistic. Finally, BFBMAiitc has the advantage that it is easily applied to 
large numbers of tissues; neither of the other methods scales well, either computationally or statistically, 
with the number of tissues, because the number of terms in the sum in equation ([s]) is 2"^ — 1. 

When there is strong evidence against Hq, the Bayes Factors can also be used to assess which al- 
ternative configurations 7 are consistent with the data. Specifically the posterior probability on each 
configuration is: 

71 BF 

P(true configuration is 7 I data, Hq false) = =^ — ^^^i (4) 

and the posterior probability that the SNP is an cQTL in tissue s is obtained by summing the probabilities 
over configurations in which 7s = 1: 

P(eQTL in tissue s \ data, Hq false) = P(true configuration is 7 | data, Hq false). (5) 

7:73=1 
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The second of these is particularly helpful when the data are informative for an eQTL in tissue s, but 
ambiguous in other tissues: in such a case the probability ([S]) will be close to 1, even though the "true" 
configuration will be uncertain (so none of the probabilities Q will be close to 1). Because both Q 
and ([5| are sensitive to choice of hyper-parameters, we compute them using BFg^^ (where the hyper- 
parameters are estimated from the data). 

Further details of methods used are provided in the Methods section. 



Simulations 

Povirer to detect eQTLs 

We begin by comparing the ability of different methods to reject the global null hypothesis Hq; i.e. to 
detect eQTLs that occur in any tissues. We expect that a tissue-by-tissue analysis, which analyzes each 
tissue separately, will perform well for detecting eQTLs that are present in a single tissue. Conversely, 
we expect joint analysis of all tissues to perform well for detecting eQTLs that are present in all tissues. 
Our Bayesian model averaging (BMA) approach attempts, by averaging over different possible eQTL 
configurations, to combine the advantages of both types of analysis, and thus aims to perform well across 
all scenarios. 

To assess this we performed a series of simulations, with five tissues measured in 100 individuals (and 
no intra-individual correlations). Each simulation consisted of 2,000 gene-SNP pairs (i.e. one candidate 
SNP per gene), half of which were "null" (i.e. the SNP was not an eQTL in any tissue), and the other 
half following an alternative hypothesis where the SNP was an eQTL in exactly k tissues, with k varying 
from 1 to 5. Thus, for example, the simulations with k = 1 assess power to detect eQTLs that are active 
in just one tissue, whereas the simulations with k — 5 assess power to detect eQTLs that are active in 
all five tissues. When simulating eQTLs that are active in multiple tissues we assumed their effects to 
be similar, but not identical, across tissues (see Methods). We applied four analysis methods to these 
data: 1) BFbmAj being our Bayesian Model Averaging approach with default weights described above; 
2) BFeMAiitc, being the computationally-scalable version of BMA described above; 3) ANOVA/linear 



regression (ANOVA/LR) (c.f. 15 and see Methods), which jointly analyzes all tissues in a regression 
model, and compares the general alternative model (which allows a different genetic effect in each tissue) 
with the null model (no effect in all tissues); and 4) a "tissue- by-tissue" analysis (c.f. [l2]), where we 
use linear regression to test for an eQTL separately in each tissue, and take the minimum p-value across 
tissues as a test statistic. For simplicity we defer consideration of the more sophisticated of our approaches, 
BFbmAi to slightly more complex simulations described later. Each of these methods produces a test 
statistic for each SNP-gene pair, testing the global null hypothesis Hq- For each test statistic, we found 
the threshold that yielded a False Discovery Rate of 0.05 (based on the known null/alternative status of 
each SNP-gene pair), and assessed the effectiveness of each method by the number of discoveries it made 
at that FDR. 

The results of these comparisons are shown in Figure [T]'V. As expected, for eQTLs that occur in just 
one tissue, the tissue-by-tissue analysis performs best. However, it is only slightly more effective than 
the joint analysis approaches in this setting. Conversely, the joint analysis approaches outperform the 
tissue-by-tissue analysis for eQTLs that occur in more than one tissue, with the gains becoming larger 
as the number of tissues sharing the eQTL increases. The BMA analyses generally perform similarly to 
one another, and outperform ANOVA/LR. This is presumably because our simulations involved eQTLs 
that have similar effects in each tissue, and our prior distribution p(l3\j,9) explicitly up-weights eQTLs 
with this feature. 

This first set of simulations assumed error variances to be equal among tissues. This assumption is 
made by ANOVA/LR, but not by the other methods, and is likely often to be violated in practice. To 
assess the effects of this we repeated the simulations, but with error variances differing among tissues. 
The results (Figure [T|3) confirm that, relative to other methods, ANOVA/LR performs less well when 
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Figure 1. The BFbma joint analysis has more power across a range of alternatives. A. Five 
tissues are simulated, each with the error variance equal to 1. B. Five tissues are simulated, with error 
variances being 1, 1.5 or 2. C. Twenty tissues are simulated, each with the error variance equal to 1. 



error variances vary among tissues. 

To assess performance in larger numbers of tissues we repeated the simulations above, but with 20 
tissues (so k = 1, . . • , 20). For this many tissues computing BFbma involves averaging over all 2^^ > 10^ 
possible eQTL configurations, which is computationally inconvenient, so we omitted BFbma from this 
comparison. The results (Figure [TjC) show that BFBMAUto performs similarly to the tissuc-by-tissue 
analysis for eQTLs that occur in just one or two tissues, and outperforms it substantially for eQTLs 
occurring in many tissues. As expected, ANOVA/LR outperforms tissue-by-tissue analysis for eQTLs 
occurring in many tissues, but is noticeably less effective for eQTLs occurring in only one tissue, and 
performs consistently less well than BFBMAiitc- 

In summary, these simulations illustrate the benefits of Bayesian Model Averaging as a general strategy 
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for producing powerful test statistics: by explicitly averaging over a range of alternative models, these 
test statistics are able to achieve good power to detect a wide range of different types of eQTL. 

Identifying eQTLs in particular tissues: borrowing information among tissues 

Next we consider the benefits of jointly analyzing multiple tissues, even when the main goal is to identify 
eQTLs in a particular tissue of interest. For intuition, suppose for a moment that every eQTL is shared 
among all tissues. Then, from the simulation results above, we know that a joint analysis will identify 
more eQTLs overall, and hence more eQTLs in the tissue of interest. Of course, not all eQTLs are shared 
among all tissues, but some are, and some tissues may share eQTLs more than others. To allow for 
this, our hierarchical model attempts to infer the extent of such sharing (by estimating the configuration 
weights 7]^), and exploits any sharing that does occur to increase power to detect eQTLs in each tissue. 
By estimating sharing from the data, rather than assuming that all tissues share equally with one another 
(as do the simpler test statistics BFbma and BPsMAUte used above), we expect BF™^ to make more 
effective use of sharing in the data to further improve power to identify eQTLs. 

To illustrate this, we simulated eQTL data for five tissues. Some eQTLs were shared by all tissues, 
some were specific to each tissue, some were shared by Tissues 1 and 2 only, and some were shared by 
Tissues 3, 4 and 5. To show how the benefits of sharing can change with sample size, we simulated 60 
samples for Tissue 1, and 100 samples for the others. This mimics a setting where Tissue 1 is harder to 
obtain than the other tissues, with Tissue 2 being the best proxy for Tissue 1. 

We applied our Bayesian methods and a tissue-by-tissue analysis to these data, and assessed their 
ability to identify eQTLs in each tissue. For the tissue-by-tissue analysis the test statistic in each tissue 
is simply the linear regression p-value in that tissue. For our Bayesian methods, the test statistic in each 
tissue is the posterior probability of the SNP being an active eQTL in that tissue ([s]). Note that this 
posterior probability is computed from joint analysis of all tissues, and takes account of sharing of eQTLs 
among tissues. For example, consider a SNP showing modest association with expression in Tissue 1. If 
this SNP also shows strong association in the other tissues, then it will be assigned a higher probability 
of being an active eQTL in Tissue 1 than it would if it showed no association in the other tissues. For 
each method, separately in each tissue, we identified the threshold of the test statistic value that yields 
a FDR of 0.05 in that tissue, based on the true active/inactive status of each SNP in that tissue (known 
since this is simulated data). (A SNP that is an eQTL in some tissues but not others counts as a "false 
discovery" if it is called as an eQTL in a tissue where it is inactive.) For the Bayesian methods we 
obtained results both using "default" weights on configurations (BFbma), and using weights estimated 
from the data by the hierarchical model (BPg^A.)- The latter is expected to be more effective as it should 
learn, for example, that Tissue 1 shares more eQTLs with Tissue 2 than with other tissues. 



The results (Figure^ show that, for all tissues, our joint analyses outperform the tissue-by-tissue 

analysis. Further, BFgf^/;^ outperforms BFbma, demonstrating the benefits of learning patterns of sharing 
from the data. Finally, the gain of BF™^ is greater for Tissue 1 than for Tissue 2, illustrating that 
benefits of sharing information are greater for tissues with small sample sizes. 

Furthermore, using the hierarchical model which pools all genes together, we can estimate the config- 
uration proportions, i.e. ry^. In the setting described above, we simulated one thousand eQTLs in each 
of 8 different configurations, as well as one thousand genes with no eQTLs. Averaged over replicates, 
the proportions are estimated to be in [0.124 — 0.127] for each of the 8 active configurations (negligible 
differences between replicates). These estimates are fairly accurate knowing that the true proportion is 
1/8 = 0.125 for each configuration. 

Analysis of eQTL data in three cell types from Dimas et al. 

We now analyze data from [12| , consisting of gene expression levels measured in fibroblasts, LCLs and 
T-cells from 75 unrelated individuals genotyped at approximately 400,000 SNPs. The data were pre- 
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Figure 2. The BF^^^ joint analysis efficiently borrows information across genes. Five 
tissues are simulated. Some cQTLs were shared by all tissues, some were specific to each tissue, and, as 
depicted by the cladogram, some were shared by Tissues 1 and 2 only, while others were shared by 
Tissues 3, 4 and 5. Each tissue has 100 samples, except tissue 1 which has only 60. 



processed similarly to the original publication, as described in the Methods section. Throughout we 
focus on testing SNPs that lie within 1Mb of the transcription start site of each gene (the "cis candidate 
region"), and on a subset of 5012 genes robustly expressed in all three cell-types. 

Gain in power from joint analysis 

First we assess the gain in power from mapping eQTLs in all three cell types jointly, using BFbma, 



compared with a "tissue-by-tissue" analysis similar to that in 12 . For each method we compute a test 



statistic for each gene, combining information across SNPs, to assess the overall support for any eQTL in 
that gene in any tissue. For our Bayesian approach the test statistic is the average value of BFbma over 
all SNPs in the cis candidate region; for the tissue-by-tissue analysis the test statistic is the minimum p- 
value from linear regressions performed separately in each tissue for each SNP in the cis candidate region. 
We translate each of these test statistics into a p-value for each gene by comparing observed values with 



9 



simulated values obtained under Hq (by permutation of individual labels). Finally we computed, for 
each method, the number of genes identified as having an eQTL at an FDR of 0.05 (using the qvalue 
package [27]). 
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Figure 3. The BFbma joint analysis is more powerful on the data set from Dimas et al. A 

and B. Histograms of gene p-values obtained by the tissue-by-tissue analysis and the joint analysis. C. 
Numbers of eQTLs called by both methods or either one of them. 



Joint mapping, via BFbma, substantially increased power to identify eQTLs compared with tissue- 
by-tissue analysis. For example, BFbma identified 1022 eQTLs at FDR=0.05, which is 63% more than 
the 627 eQTLs identified by the tissue-by-tissue analysis at the same FDR (Figure |3] A and B). Further, 
the vast majority of eQTLs identified by the tissue-by-tissue analysis (94%) are also detected by BFbma 
(Figure [Sp). 

In many cases, the eQTLs detected by BFbma but not by the tissue-by-tissue analysis have modest 
effects that are consistent across tissues. Because their effects are modest in each tissue, they fail to 
reach the threshold for statistical significance in any single tissue, and so the tissue-by-tissue analysis 
misses them. But because their effects are consistent across tissues, the joint analysis is able to detect 
them. Figure |4] shows an example illustrating this (gene ASCCl, Ensembl id ENSG00000138303, with 
SNP rsl678614). The PC-corrected phenotypes already indicate that this gene-SNP pair looks like a 
consistent eQTL (Figure |4|\) , and its effect size estimates are highly concordant across tissues (Figure 
[4]B). However, as indicated by the g- values on the forest plot, this eQTL is not called by the tissue-by- 
tissue analysis in any tissue (all the g- values exceed .14). In contrast, the joint analysis pools information 
across tissues to conclude that there is strong evidence for an eQTL [q — 0.001), and that it is likely an 
eQTL in all three tissues (probability 1 assigned to the consistent configuration 7 = [111], conditional on 
it being an eQTL). 



Many eQTLs are consistent among tissues 

The original analyses of these data concluded that 69% to 80% of eQTLs operate in a cell-type specific 
manner ( [l2]). These results were obtained by mapping eQTLs separately in each tissue, and then exam- 
ining which of the eQTLs identified in each tissue also showed some signal (e.g. at a relaxed significance 
threshold of p = 0.05), in another tissue. However, as noted by 14 , due to incomplete power, eQTLs that 
are actually shared between tissues may appear "tissue-specific" in this type of analysis. Our hierarchical 
model has the potential to help overcome this difficulty by estimating the proportion of eQTLs that 
follow each configuration type as a parameter of the model, combining information across all genes, and 
without setting specific significance thresholds (thus sidestepping the problems of incomplete power). 
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Figure 4. Example of an eQTL with weak, yet consistent effects. A. Boxplots of tlie 
PC-corrected expression levels from gene ASCCl (Ensembl id ENSG00000138303) in all three cell 
types, color-coded by genotype class at SNP rsl678614. B. Forest plot of estimated standardized effect 
sizes of this eQTL. Note that none of the q-values from the tissue-by-tissue analysis are significant at 
FDR=G.05. 



Applying the hierarchical model to these data produced an estimate of just 8% of eQTLs being specific 
to a single tissue, with an estimated 88% of eQTLs being common to all three tissues (95% CI = 84%- 
93%; Tablejl]). Among eQTLs shared between just two tissues, many more are shared between LCLs and 
T-cells, than between these cell types and fibroblasts. This is consistent with results from [l2], and not 
unexpected since LCLs and T-cells are more similar to one another than to fibroblasts. 

Table 1. Inference of the proportion of tissue specificity 



Configuration 
F-L-T 
L-T 
F-L 
F-T 
F 
L 
T 



Hierarchical 
0.882 [0.840, 
0.051 [0.025, 
0.005 [0.000, 
0.002 [0.000, 
0.033 [0.014, 
0.015 [0.000, 
0.011 [0.000, 



model 
0.925] 
0.085] 
0.018] 
0.011] 
0.065] 
0.039] 
0.033] 



Tissue-by-tissue 
0.187 
0.080 
0.050 
0.047 
0.246 
0.165 
0.224 



The results for the hierarchical model were obtained with the multivariate Bayes Factors allowing 
correlated residuals and the EM algorithm. The results for the tissue-by-tissue analysis were obtained 
by calling eQTLs at an FDR of 0.05 after performing permutations in each tissue separately, and 
calculating the overlaps among tissues. 



We obtained qualitatively similar patterns when we varied some of the assumptions in the hierarchical 
model - specifically, whether or not we allow for intra-individual correlations, whether or not we assume 
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at most one eQTL per gene, whether or not we remove PCs to account for confounders, and whether or 
not we analyze ah genes or only those genes robustly expressed in all tissues (Supplementary text SI). 
Nonetheless, we caution against putting too much weight on any particular number to quantify tissue 
specificity, not least because the definition of a tissue-specific eQTL is somewhat delicate: for example, it 
is unclear how to classify a SNP that is very strong eQTL in one tissue, and much weaker in the others. 
Further, our estimates necessarily reflect patterns of sharing only for moderately strong cQTLs, strong 
enough to be detected in the modest sample sizes available here: patterns of sharing could be different 
among weaker eQTLs. Nonetheless, these results do suggest that there is substantial sharing of eQTLs 
among these three tissue types, considerably higher than the original analysis suggested. 

To illustrate the potential pitfalls of investigating heterogeneity in a tissue-by-tissue analyses, we also 
ran a tissue-by-tissue analysis on these data. Specifically, we called eQTLs separately in each tissue 
(at an FDR of 0.05), and then examined the overlap in the genes identified in each tissue. Using this 
procedure, in strong contrast with results from the joint analysis, 65% of eQTLs are called in only one 
tissue, with fewer than 15% called in all three tissues (Table [T]). Qualitatively similar results are obtained 
for different FDR thresholds. However, these results cannot be taken as reliable indications of tissue 
specificity, because the procedure fails to take account of incomplete power to detect eQTLs at any given 
threshold, and therefore tends to over-estimate tissue specificity. Figure [5] shows an eQTL that illustrates 
this behavior (gene CHPTl, Ensembl id ENSG00000111666, with SNP rsl0860794). Visual examination 
of the expression levels in each genotype class (Figure [sJfV), suggest that this SNP is an eQTL in all three 
tissues, with similar effects in each tissue (Figure [5|3). This is supported by the joint analysis, which 
shows strong evidence for an eQTL q = 0.001, and assigns probability effectively 1 to the consistent 
configuration 7 — [111]. However, as shown by the g- values , at an FDR of 0.05, the tissue-by-tissue 
analysis calls the eQTL only in fibroblasts. 
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Figure 5. Example of an eQTL wrongly called as tissue-specific by the tissue-by-tissue 
analysis. A. Boxplots of the PC-corrected expression levels from gene CHPTl (Ensembl id 
ENSG00000111666) in aU three cell types, color-coded by genotype class at SNP rsl0860794. B. Forest 
plot of estimated standardized effect sizes of this eQTL. Note that, from the g- values of the 
tissue-by-tissue analysis, the eQTL is significant at FDR=0.05 only in fibroblasts. 
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Given the disagreement between the results from our novel framework and the original analyses of 
these data, we checked the plausibility of our results by applying a previously-used method for examining 
pairs of tissues to these data ( (13j). This analysis takes the best eQTL in each gene identified in 
one tissue, and then estimates the proportion of these (tti) that are also eQTLs in a second tissue (by 



applying Storey's method 27 to their nominal p-values in that second tissue, uncorrected for multiple 
comparisons). Unlike the tissue-by-tissue analysis above, this approach avoids thresholding of p- values , 
and makes some allowance for incomplete power. However, unlike our framework, this approach can only 
be applied to compare pairs of tissues. Applying this approach to each pair yielded a mean estimate 
of TTi « 88% (range 77% to 94%), broadly consistent with our qualitative conclusion that a substantial 
proportion of cQTLs are shared among tissues. 



Discussion 

In this work, we have presented a statistical framework for analyzing and identifying eQTLs, combining 
data from multiple tissues. Our approach considers a range of alternative models, one for each possible 
configuration of eQTL sharing among tissues. We compute Bayes Factors that quantify the support in 
the data for each possible configuration, and these are used both to develop powerful test statistics for 
detecting genes that have an eQTL in at least one tissue (by Bayesian model averaging across configu- 
rations), and to identify the tissue(s) in which these eQTLs are active (by comparing the Bayes factors 
for different configurations against one another) . Our framework allows for heterogeneity of eQTL effects 
among tissues in which the eQTL is active, for different variances of gene expression measurements in 
each tissue, and for intra-individual correlations that may exist due to samples being obtained from the 
same individuals. For eQTL detection, our framework provides consistent, and sometimes substantial, 
gains in power compared to a tissue-by-tissue analysis and ANOVA or simple linear regression. Con- 
cerning the tissue specificity of eQTLs, our framework efficiently borrows information across genes to 
estimate configuration proportions, and then uses these estimates to assess the evidence for each possible 
configuration. When re-analyzing the gene expression levels in three cell types from 75 individuals ( ^ 12] ) , 
we found that there appears to be a substantial amount of sharing of cQTLs among tissues, substantially 
more than suggested by the original analysis. 

In the next few years, we expect that expression data will be available on large numbers of diverse 
tissue types in sufficient sample sizes to allow eQTLs to be mapped effectively (for example, the NIH 
GTEx project aims to collect such data). The methods presented here represent a substantive step 
towards improved analyses that fully exploit the richness of these kinds of data. However, we also see 
several directions for potential extensions and improvements. First, our current framework can only 
partially deal with the challenges of large numbers of tissues. Specifically, because with S tissues, there 
are 2^^ possible configurations of eQTL sharing among tissues, some of our current methods, which 
consider all possible configurations, will become impractical for moderate S (speculatively, above about 
10, perhaps). Our test statistic BFeMAUto partially addresses this problem, by allowing for heterogeneity 
while averaging over only 5-1-1 configurations, which is practical for very large S. Our simulation 
results suggest that BFeMAUto is a powerful test statistic for identifying SNPs that are an eQTL in at 
least one tissue. However our preferred approach for identifying which tissues such SNPs are active in 
involves a hierarchical model that estimates the frequency of different patterns of sharing from the data, 
and this hierarchical model scales poorly with S. In particular, having a separate parameter for each 
possible configuration is unattractive (both statistically and computationally) for large S, and alternative 
approaches will likely be required. There are several possible ways forward here: for example, one would 
be to reduce the number of distinct configurations by clustering "similar" configurations together; another 
would be to focus less on the discrete configurations, and instead to focus on modeling heterogeneity in 
effect sizes in a continuous way - perhaps using a mixtures of multivariate normal distributions with more 
complex covariance structures than we allow here. We expect this to remain an area of active research 
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in the coming years, especially since these types of issues will likely arise in many genomics applications 
involving multiple cell types, and not only in eQTL mapping. 

Another important issue to address is that most future expression data sets will likely be collected by 
RNA-seq, which provides count data that are not normally distributed. Previous eQTL analyses of RNA- 
seq (e.g. f28l) have nonetheless performed eQTL mapping using a normal model, by first transforming 
(normalized) count data at each gene to the quantiles of a standard normal distribution. Although this 
approach would not be attractive in experiments with small sample sizes, with the moderate to large 
sample sizes typically used in eQTL mapping experiments this approach works well. As a first step, this 
approach could also be used to apply our methods to count data. However, ultimately it would seem 
preferable to replace the normal model with a model that is better adapted to count-based data, perhaps 
a quasi-Poisson generalized linear model ( [29]); Bayes Factors under these models could be approximated 
using Laplace approximations, similar to the approximations used here for the normal model 19 . The 
quasi-Poisson model has the advantage over the normal transformation approach that it preserves the 
fact that there is more information about eQTL effects in tissues where a gene is high expressed than 
in tissues where it is low expressed. This information is lost by normal transformation. In our analyses 
here we addressed this by analyzing only genes that were robustly expressed in all tissues, but this is 
sub-optimal, and will become increasingly unattractive as the number of tissues grows. 



Methods 

Materials and Methods 

Software implementing our methods are available on the website http : / / stephenslab . uchicago . edu/ software . html. 

Bayesian Methods for Mapping Multiple-tissue eQTLs 
Models for Multiple-tissue eQTLs 

For each tissue, we model the potential genetic association between a target SNP and the expression 
levels of a target gene by the simple linear regression model ([!]) . In vector form, this model is represented 
by 

y, = M.l + PsQs + e., e, - U{0, all), (6) 

where s indexes one of the S tissue types examined and the vectors VstOs ^iid denote the expression 
levels, the genotypes of the samples and the residual errors respectively for the s*^ tissue type. The 
intercept term, /is, and the residual error variance, are allowed to vary with tissue type. The regression 



coefficient /3s denotes the effect of the eQTL in tissue s, but we follow 19 24 in using the (unitless) 
standardized regression coefficient := /Ss/cts, as the main measure of effect size. As a result, inference 
is invariant to scale transformations of the response variables (j/s) within each tissue. 

When the tissue samples are taken from the same individuals we allow that the observations on the 
same individual may be correlated with one another. Specifically, let E := (ei • • • Cs) denote the N x S 
matrix of residual errors, the we assume it to follow a matrix- variate normal (MN) distribution, i.e., 

£; - MN(0,/,S). (7) 

That is, the vectors {en, . . . , esi) are independent and identically distributed as A/'(0, S). The (unknown) 
S X S covariance matrix E quantifies the correlations between the S tissues; it can vary from gene to 
gene and is estimated from the data (see below). [When the tissue samples are collected from different 



individuals then we assume their error terms are independent; methods for this case are given in 19 .] 



14 



Prior on effect sizes 

A key component of our Bayesian model is the distribution p{b\^,9), where 6 denotes hyper-parameters 
that are to be specified or estimated from the data. (In the main text we used p{/3\"/, 9) to simphfy 
exposition, but we actually work with the standardized effects b.) Of course, if 7s = then &s = by 
definition. So it remains to specify the distribution of the remaining bs values for which 7^ = 1. 



We use the distribution from 19 (see also 17 18 ), which provides a flexible way to model the 



heterogeneity of genetic effects of an eQTL in multiple tissues. Specifically, 19 consider a distribution 
p{b\<f), uj, 7), with two hyper-parameters, 0, w, in which the non-zero effects are normally distributed about 
some mean &, which itself is normally distributed: 

6s|6,7. = l^AA(6,</.2), (8) 

and 

6-A/'(0,tj2). (9) 

Note that (fi^ + uP controls the variance (and hence the expected absolute size) of 6s, and (fp' /{(fp' -I- w^) 
controls the heterogeneity (indeed, oj^/ {(fp -f w^) is the correlation of hg, bs' for different subgroups s ^ s'). 
If (fp ^ then this model corresponds to the "fixed effects" model in which the effects in all subgroups 
are equal (e.g. [20]). 

To allow for different levels of effect size and heterogeneity, jl9] use a fixed grid of values {((/), , w,) : 
i = 1, . . . , L}, with the ith grid point having weight Wi. Thus 

p(6|7,^?) =^i(;,p(6|0„w„7). (10) 

i 

In all our applications here we consider the grid of values fixed, and treat the weights Wi, . . . ,wl as 
hyper-parameters (so 9 = {wi, . . . , wl)), which can be either fixed or estimated. 



Choice of grid for ((/>, a;) 

We define a grid of values for {4>,uj) by specifying a set E of values for the average effect size, uj^ + (jP, 
and a set H of values for the heterogeneity (jP /{(fP -f w^), and then taking the grid to be all L = \E\x \H\ 
possible combinations of values. For all methods we use E — {0.1^, 0.2^, 0.4^, 0.8^, 1.6^}, which is designed 
to span a wide range of eQTL effect sizes. For BFbma and BFg^^ we allow for only a limited range of 
heterogeneity: H = {0,0.25}. In this way we assume that when the eQTL is present in multiple tissues, 
it has a similar (but not necessarily identical) effect in each tissue. For BFeMAUte we allow a much wider 
range of heterogeneity: H = {0,0.25,0.5,0.75,1}. The rationale here is that the large heterogeneity 
values will help capture eQTLs that are present in only a subset of tissues, a feature that is not otherwise 
captured by BFeMAiitc as it averages over a small number of configurations. 



Choice of weights w and 

Let I7I := 7s denote the nmnber of elements of 7 that are equal to 1 (i.e. the number of tissues in 
which the eQTL is active in configuration 7). 

For BFbma we fix the weights so that they put weight 1/5 on all S possible values for I7I, and, 
conditional on I7I, put equal weight on all (|^|) configurations with that value for |7|. Thus rj-^ = 

(1/S')(|^l)" . In addition we fix the grid weights w to be equal on all grid values. 

For BFeMAiitc we put non-zero weights ry on only the consistent configuration (I7I = S) and configu- 
rations with an eQTL in a single tissue (I7I = 1). We set 77 so that it puts weight 0.5 on each of I7I 1 
and I7I = S. Conditional on I7I = 1 we assume all S possibilities are equally likely. Thus rj^ — 0.5 if 
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7 = [111 ... 1] and 0.5/S if I7I = 1. Again, we fix the grid weights w to be equal on aU grid values (but 
with the larger grid for heterogeneity described above). 

For BFg^A estimate the weights w,r) from the data using a hierarchical model to combine infor- 
mation across genes, as described below. 

Bayes Factor Computation 

To complete model specification, we use (limiting, diffuse) prior distributions for the nuisance parameters 
fis and E, as in ^30J . Under these priors we can compute the Bayes Factor BF-y in ([2]) using 

M 

BF^ = ^u;,BF^(0„a.,) (11) 

where BF^(0j,u;j) is given by 

^ p{Y\G,(j},uj,j) ^ J p{Y\G,^i,b,j:)p{^i,j:)p{b\j,(p,uj)dbd^dJ: 

piY\G,Ho) Jp{Y\G,^i,b^O,i:)p{(i,i:)dfid^ ^ ' 

where Y and G denote the collection of expression levels and genotypes for a target gene-SNP pair across 
all tissue types respectively. We use analytic approximations for these Bayes Factors based on Laplace 



approximation, given in 19 30 . In particular, we use the approximation which in functional forms is 
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which combines information across genes, 



connected to Frequentist's score statistic. 

Bayesian Hierarchical Model 

For BFg^j(^ we use a hierarchical model, similar to 
to estimate the grid weights w^s and configuration weights jy's. Following both 21 22 we make the 
simplifying assumption that each gene has at most one eQTL (which may be active in multiple tissues), 
and that each SNP is equally likely to be the eQTL. Let rm. be the number of SNPs in the cis-region for 
gene k. Then, if BF!^'"(0, oj) denotes the Bayes Factor il2n computed for SNP v in gene k, the "overall 
Bayes Factor" measuring the evidence for an eQTL in gene k, BF*^, is obtained by averaging over the 
possible eQTL SNPs, the possible configurations 7, and the grid of values for 4>,uj, weighting by their 
probabilities: 

T^T^fr/ N p(data at gene fclgene contains eQTL) ,^ , xv-^v^v^ ^-r^k ,w : n 

= p(data at gene fclgene contains no eQTL) ^ (V-.) E E E '^.-^BF-(4, ^.). (13) 

Furthermore, if we let ttq denote the probability that each gene follows the null (i.e. contains no eQTL) 
then the likelihood for gene /c, as a function of 7To,ri,w, is given by 

Lk{T^o,r},w) = (1 — 7ro)p(data at gene fcjgene contains eQTL) + 7rop(data at gene /c|gene contains no eQTL) 

(14) 

cx (1 - 7ro)BF'= + vro (15) 
The overall likelihood for our hierarchical model is obtained by multiplying these likelihoods across genes: 

L{-K^,ri,w) ^^Lk{TiQ,'n,w). (16) 

k 

Note that although the expression levels for different genes are not independent, because the SNPs being 
tested in different genes are mostly independent this independence assumption for the likelihoods across 
genes is a reasonable starting point. We have developed both an EM algorithm to estimate the parameters 
(7ro,Tj,i(;) by maximum likelihood (see Supplementary information). 
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Relaxation of "one cis-eQTL per gene" Assumption 

To relax the "one cis-eQTL per gene" assumption we adopt the following procedure. First we compute 
the posterior probability of each SNP being the sole eQTL for each gene (i.e. only allowing one cis-eQTL 
per gene) with a set of default parameters, and use these to identify the top SNP for each gene (i.e. the 
one with the largest posterior probability of being the eQTL). 

For each gene, separately in each tissue, we compute the residuals of its expression level after regressing 
out the effect of the top SNP. If these residuals are strongly associated with a SNP then this is evidence 
for that SNP being a second independent eQTL for that gene. Therefore, to allow for the potential for 
more than one eQTL per gene we treat these residuals as defining a second set of "artificial" expression 
data for each gene and each tissue, and fit the hierarchical model using both the original and the artificial 
expression data. 



Simulation procedures 

For our simulations, when simulating SNP-gene pairs, the genotypes at each SNP in each individual were 
simulated as Binomial(2,0.3): that is, with minor allele frequency 30% and assuming Hardy- Weinberg 
equilibrium. Phenotypes with eQTLs were simulated, with effect size based on an expected proportion 
of variance explained (PVE) of 20%; (see supplementary text SI). For Figures ^ and [l^) the error 
variances (one per tissue) were all equal to 1. For Figure [Ip the error variances were randomly drawn 
from {1, 1.5, 2}, all equally likely. 



The ANOVA/LR method 

The ANOVA/LR method uses the same linear model as our Bayesian methods ([T]), except that the 
residual errors <Js are assumed to be equal across tissues s. Within this model we tested the global null 
hypothesis (/3s = for all s) using an F test comparing the null model with the unconstrained alternative 
(/3s unconstrained). 



Preprocessing of the data set from Dimas et al. 

The phenotypes from Dimas et al. ( [l2]) were retrieved from the Gene Expression Omnibus (GSE17080). 
We mapped the 22,651 non-redundant probes to the hgl9 human genome reference sequence (only the 
autosomes) using BWA ( 31 1), kept 19,965 probes mapping uniquely with at most one mismatch, and 
removed the probes overlapping several genes from Ensembl. This gave us 12,046 genes overlapped by 
16,453 probes. For genes overlapped by multiple probes, we chose a single probe at random. In our 
analyses we considered only genes that were robustly expressed in all tissues. A gene was considered 
robustly expressed in a given tissue if its mean expression level across individuals in this tissue was larger 
than or equal to the median expression level of all genes across all individuals in this tissue. As a result, 
we focused on 5012 genes. 

Genotypes were obtained from the European Genome-phenome Archive (EGAD00000000027). We 
extracted the genotypes corresponding to the 85 individuals for which we had phenotypes and converted 
the SNP coordinates to the hgl9 reference using liftOver ( [32]). To detect outliers, we performed a 
PGA of these genotypes using individuals from the GEU, GHB, JPT and YRI populations of the HapMap 
project using EIGENSOFT ( [33j). As in the original study, we identified 10 outliers and removed them 
from all further analyses, which were therefore performed on 75 individuals. 

Gene expression measurements suffer from various confounders, many of which may be unmeasured 
( [34|), but which can be corrected for using methods such as principal components analysis (PGA). 



Following 28 , we applied PGA in each tissue separately on the 5012 x 75 matrix of expression levels 



of each gene in each individual. We sorted principal components (PGs) according to the proportion of 
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variation in the original matrix they explain, and selected PCs so that adding another PC would explain 
less than 0.0025% of the variation. As a result, this procedure identified 16 PCs in Fibroblasts, 7 in LCLs 
and 15 in T-cells. We then regressed out these PCs from the original matrix of gene expression levels, 
and used the residuals as phenotypes for all analyses. 

All methods we compared assume that the errors are distributed according to a Normal distribution. 
Before analysis we therefore rank-transformed the expression levels at each gene to the quantiles of a 
standard Normal distribution ( [24]). 



Permutation procedures 

On the data set from Dimas et al., we assessed the performance of two methods, the tissue-by-tissue 
analysis and the BMA joint analysis, by comparing the number of genes identified as having at least one 
eQTL in any tissue, at a given FDR. For each method, we defined a test statistic, which was computed 
for each gene. For the tissue-by-tissue analysis, the test statistic is the minimum p-value of the linear 
regressions between the given gene and each cis SNP in each tissue (so the minimum is taken across all 
SNPs and all tissues). For the BMA joint analysis, the test statistic is the average of the Bayes Factors 
for the given gene and each cis SNP. (When applying the tissue-by-tissue analysis to test for eQTLs in 
a single tissue, the test statistic is the minimum p- value of the linear regressions between the given gene 
and each cis SNP in that tissue.) 

In each case we converted the test statistic to a p-value for each gene, testing the null hypothesis 
that the gene contains no eQTL in any tissue, by comparing the observed test statistic with the value 
of the test statistic obtained on permuted data obtained by permuting the individuals labels (using the 
same permutations in each tissue to preserve any intra-individual correlations between gene expression 
in different tissues). Specifically, let P denote the total number of permutations (we used P = 10*), 
Tg the value of the test statistic for gene g on the non-permuted data, and Tg the value of the test 
statistic on the i*''-permuted data. The p-value for gene g from the tissue-by-tissue analysis is: (1 -I- 
J^Li 1t(-)<t )/(1 + P)- For the BMA joint analysis, the p-value is: (1 + J^Li K(')^t )/(1 + P)- Note 

^ g -^^ g ^ g g 

that permutations were performed for each gene, since the null distribution of the test statistic will vary 
across genes (not least because the genes have different numbers of SNPs in their cis candidate region; 
see supplementary figure SI). 

From the p-value calculated for each gene we estimate q values [2^ using the qvalue package, and 
determine the number of genes having at least one eQTL in any tissue at an FDR of a by computing the 
number of genes with q < a. 

When performing the tissue-by-tissue analysis on a single tissue, we performed the permutations in 
each tissue separately. 
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A Computational algorithm for fitting hierarchical model 

For the hierarchical model described in the main text, our primary interest is making inference on the 
parameter set Q = (ttq,??, A). Here, we give details of an algorithm for inferring O, via maximum 
likelihood estimation based on the EM algorithm. 
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A.l Notations 

Throughout this section, we adopt the following additional notations. For gene k, we use a latent binary 
indicator Zk to denote if there is any eQTL in its cis-region for any tissue type, in particular, 

P{zk = 1) = 1 - TTo; (17) 

Wc use a latent random indicator mfe-vcctor Sfe to denote the true eQTL SNP conditional on = 1 and 
let Skp denote the p-th entry of s^. The "one cis eQTL per gene" assumption restricts can have at 
most one entry equaling 1 (with the remaining entries being 0). By this definition, 

P(sfe=O|zfe = 0) = l, (18) 

and we also make the simplifying assumption that 

P{skp = Ikfe = 1) = — . (19) 

Furthermore, for gene k and SNP p, we index all configurations and use a (2"^ — l)-dimension latent 
indicator vector Ckp to denote the actual configuration for the gene-SNP pair. In case the SNP is not an 
eQTL, 

P{ckp = 0\skp = 0) = 1. (20) 
Otherwise, we assume the jth configuration is active with prior probability 

P(cfepj l|sfep = 1) = Vj- (21) 

Joining the column vectors Ckp for all ruk SNPs, we obtain a latent (2'^ — 1) x rUk random matrix C^. 
Finally, we use the latent i- vector Wkp indicate the actual prior effect size for active tissue types for the 
pair of gene k and SNP p. The m-th entry of the indicator is denoted by Wkpm, for which we assume 
prior probability 

P{wkp = 0\skp = 0) = 1, (22) 

and 

Piwkpm = l\skp = 1) = Xm- (23) 

Joining the column vectors Wkp for all SNPs, we obtain a latent L x ruk random matrix Wk 
A. 2 Mctximum Likelihood Inference based on EM algorithm 

In the maximum likelihood framework, we treat latent variables Zk, Sk, Ck and Wk, k = 1, . . . ,G as missing 
data and apply the EM algorithm. 

For a total number of G genes, let z = {zi, . . . , Zg),t = {si, . . . , Sg),C = (Ci,...,Cg) and 
W = {Wi, . . .Wg) denote the complete collection of latent variables. Let Y = {Yi, . . .Yg) and 
G = (Gi, .... Gq) denote the complete set of observed data. Based on the hierarchical model described 
in previous section, wc can write out the complete data log-likelihood as follows, 

log p{Y,z,T,C,W\G,e) = 

^{1- Zk) log TTo + '^Zk l0g(l - TTo) 

k k 

+ ^ ZkSkp log — + ^ ZkSkpCkpj log IJj + ^ ZkSkpWkpm log Xm ^'^^^ 

t , . , 

kyp k,p,3 k,p,m 

+ ^ ZkSkpCkpjWkpm ■ BF/jpjTO + ^ logPfe. 
k,p,j^m k 
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In (24), p1 denotes the likelihood of the null model for gene k, i.e., 

pl:^p{Yk\zk^Q) (25) 

and 

o-p P{Yk\zk = 1, Skp = 1, Ckpj = 1, Wkprn = 1, Gfe, 9) 

0-rfcp,„i — [Zb] 

Pk 

is the Bayes Factor (pre-)computed for a fully specified alternative model. 

The EM algorithm searches for maximum likelihood estimate of 9, by iteratively performing an 
expectation (E) step and a maximization (M) step. 

In the E-step, for the t-th iteration, we evaluate the expectation of complete data log-likelihood (24) 
conditional on current estimate of parameter G and Y . The computation is straightforward, for 

example, 

E(zfc|rfe,Gfe,e(*)) = p{zk - i|rfc,Gfe,e(*)) 

^ p(zfc = i|9(*))-p(yfc|zfc^i,Gfc,e(*)) 

p{Yk\GkM'^) (27) 
(l-4*))BFi*) 



rg + (,1 - tTq )Bi'j. 



similarly. 



E(.,.,,|r,,G,,eW) = ° (28) 

E(zfc5fc,Cfep,z/;fc,^|r,G,e(*)) = ^ " " .y" ' (29) 

<' + (!-<' )BFi*> 



where 



(,) p{Yk\zk^l,GuM'^) 



V — r7(*)A(*)BF. 



(30) 



and 



(t) p(Yfe|zfe = l,Sfcp = l,GA;,e) 
B^ fcp - ^^0 



(31) 



In the M-step, we find a new set of estimates, 0^^"+^^, by maximizing the conditional expectation 
E (logp(y, z, T, G, W\G, Q)\Y, G, 6^). In this case, the simultaneous maximization can be performed 
analytically. In particular, 

1 9 (t) 

_ i V ^ r32) 

/ 7fcp A„ Bl't^p^,„, (t) 
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and 



\(t+l) _ '^'-^ ^0 +(1-^0 JBl-fc /n.x 



Typically, we initiate the EM algorithm by setting 0*^°^ to some random values and running iterations 
until some pre-defined convergence threshold is met (In practice, we monitor the increase of the the 
log-likelihood function between successive iterations, and stop the iterations as the increment becomes 
sufficiently small.). 

We construct profile likelihood confidence intervals for estimated parameters. For example, a (1 — a)% 
profile likelihood confidence set for ttq is built as 

{tto : logp(V|7ro, f], A, G) > logp(V|7ro, f], A, G) - ^^fi_,)}, (35) 

where ttq, ^, A are MLEs obtained from the EM algorithm. 



B Supplements for the simulations 

B.l Simulate eQTL data via the proportion of variance explained 

For a given gene-SNP pair at a time, we simulate data in S tissues according to a particular configuration. 
In a given tissue s G {1, ... ,5*} for which the SNP is an eQTL (/3s 7^ 0), let's define the proportion of 
variance in phenotype explained by the genotype: 



When working with standardized effect sizes bs = Ps/o's- 

V{Xbs) 



PVE.. 



As stated elsewhere ( [?]), we approximate the expectation of the PVE via a ratio of expectations, noted 

E[V{Xbs)] 
E[V{Xb,)] + 1 

We assume that the genotypes are drawn from a Binomial distribution with parameters 2 and /, 
the minor allele frequency, so that i?[y(X)] = 2/(1 — /). Moreover, as we assume 6s|6 ~ N(b,(p'^) and 
b ~ N{0,u'^), the marginal effect size is 6s A''(0, 0^ +w^). We can hence approximate b^ by its variance. 
Therefore: 

(0^+a;^)x2/(l-/) 
(</.2+a;2)x2/(l-/) + l 

By fixing h (eg. 20%) as well as the minor allele frequency (eg. 30%), we obtain: 



(1 - ft) X 2/(1 - /) 



Now if we fix the heterogeneity in effect sizes (eg. (jp' /{(jp +0;^) = 20%), we can deduce (jp and then a;^. 
We can hence draw b and then each 6,16. 
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Once we have them, it is straightforward to simulate the phenotype of the i*^ individual in the s*'' 
tissue: 

yis = bs(Jsg% +7V(0, erf) 

with Us being fixed at 1 for instance. 

B.2 Implement the ANOVA/LR model in R 

For each gene-SNP pair, the expression levels from all N individuals in all 5* tissues are recorded into a 
vector y of length N x S. The genotypes are appropriately repeated S times into a vector xg, and the 
tissue indicators are appropriately recorded into a vector xs. We can then use the ANOVA/LR model to 
test if there is an effect of the genotype with the following commands: 

ml <- lm(y ~ xs) 

ni2 <- lm(y ~ xs * xg) 

pval <- anovaCml, m2)[[6]][2] 



B.3 Calculate the empirical FDR from simulated eQTL data 

For a simulated data set of G gene-SNP pairs, let Zg be the test statistic of a given pair with g = 1, . . . , G. 
For the tissue-by-tissue method, we take as test statistic the minimum P-value across tissues. For the 
Bayesian method, the test statistic is the Bayes Factor. For the ANCOVA, the test statistic is the P-value 
of the genotype effect with interaction. 

All gene-SNP pairs can be classified as in the following table ( [27]): 

Called eQTL Not called Total 

True null F Go - F Go 

True eQTL T Gi - T Gi 

Total S G - S G 

As we simulate data, we know which pairs are true eQTLs. By fixing the empirical false discovery 
rate (FDRf.) at 5%, we can find the corresponding cutoff c on the test statistics, and from there calculate 
the true positive rate (TPR) at this cutoff: 

TPR{c) = T(c)/Gi with c such that FDRe{c) = i^(c)/5(c) = 0.05. 

The following algorithm describes how to iteratively find the cutoff c corresponding to the 5% empirical 
FDR: 
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Data: test statistics zi,...,zg 

if P-values then 
1^ sort in increasing order: Z(^i) < ■ ■ ■ < Z(^g) 

else if Bayes Factors then 
|_ sort in decreasing order: ^(i) > . . . > -2(g) 

foreach gene-SNP pair 5 1 to G do 

s -^number of called eQTLs at this cutoff c 
/ ^number of false positives among them 
fdr -(- f/s 
if fdr > 5% then 

t <— number of true positives among the called eQTLs 

tpr ^ t/s 

exit 



Between different methods, the empirical FDRs will always be 5% (or slightly higher) but the TPRs 
and FPRs will be different, which allows us to compare the performance of the methods. 

C Supplements for the analysis of the Dimas et al. data set 

C.l Hierarchical model fed with Bayes Factors from residuals 

First we computed the Bayes Factors for each combination of grid values and configurations, one gene- 
SNP pair at a time. Second, for each gene, we regressed out the effect of its best SNP, and we recomputed 
the Bayes Factors for the remaining SNPs using the residuals as phenotypes. Third, we launched the 
hierarchical model with only the Bayes Factors obtained from the residuals. 

If the "at most one eQTL per gene" assumption is reasonable for this data set, we would expect the 
estimated ttq to be very high (meaning that the vast majority of genes have no eQTL), the lowest grid 
value to have the highest probability (meaning that the effect sizes are very small), and the credible 
intervals for the configurations to be very large (corresponding to high uncertainty). 

This is indeed what we observe: 

ttq: 0.963 [0.946, LOOO] 



Grid value {(fp'^uj'^) 


Posterior mean 


95% credible interval 


(0.01, 0.01) 


0.930 


[0.543, 1.000] 


(0.01, 0.04) 


0.070 


[0.000, 0.459] 


(0.01, 0.16) 


0.000 


[0.000, 0.107] 


(0.01, 0.64) 


0.000 


[0.000, 0.036] 


(0.01, 2.56) 


0.000 


[0.000, 0.019] 
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Configuration 


Posterior mean 


95% credible interval 


100 


0.316 


[0.000, 1.000] 


010 


0.330 


[0.000, 1.000] 


001 


0.027 


[0.000, 0.535] 


110 


0.250 


[0.000, 1.000] 


101 


0.020 


[0.000, 0.442] 


Oil 


0.029 


[0.000, 0.535] 


111 


0.028 


[0.000, 0.400] 



C.2 



Configuration proportions from all genes without removing expression 
PCs 



Similarly to what was done in the first analysis of this data set ( 12 ), we also analyzed the data set 
comprising all 12,046 genes, i.e. without pre-selecting genes robustly expressed in all three tissues, and 
without removing expression PCs. Here are the configuration proportions estimated by the EM algorithm: 



Configuration 


Hierarchical model 


F-L-T 


0.793 [0.722, 0.878] 


L-T 


0.071 [0.030, 0.129] 


F-L 


0.000 [0.000, 0.015] 


F-T 


0.000 [0.000, 0.015] 


F 


0.052 [0.000, 0.109] 


L 


0.058 [0.016, 0.115] 


T 


0.025 [0.000, 0.068] 



They are thus qualitatively similar to those obtained on the subset of genes robustly expressed in all 
three tissues and after having removed PCs (table 1 of the main text). 
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Figure SI. Distribution of the number of SNPs in the cis region of each gene. 



